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This paper is concerned with classes of models of stochastic reaction dynamics with time-scales separation. 
We demonstrate that the existence of the time-scale separation naturally leads to the application of the averaging 
principle and elimination of degrees of freedom via the renormalization of transition rates of slow reactions. 
The method suggested in this work is more general than other approaches presented previously: it is not limited 
to a particular type of stochastic processes and can be applied to different types of processes describing fast 
dynamics, and also provides crossover to the case when separation of time scales is not well pronounced. We 
derive a family of exact fluctuation-dissipation relations which establish the connection between effective rates 
and the statistics of the reaction events in fast reaction channels. An illustration of the technique is provided. 
Examples show that renormalized transition rates exhibit in general non-exponential relaxation behavior with a 
broad range of possible scenarios. 
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L INTRODUCTION 

Chemical reaction networks are systems of molecular species of different types interacting with 
each other by means of multiple reactions |Ql. In classical chemical systems, the volume of the 
reactor and population numbers of species of each types are usually large giving the accurate de- 
scription of the system in terms of the concentrations. Reactors with complex chemistry give rise 
to complicated systems of nonlinear equations for the concentrations of chemical species that do 
not lend themselves to analytic solution. Dynamics of these quantities can be modeled via sets 
of ordinary differential equations (ODEs) which are powerful tools for predicting the dynamical 
behavior of macroscopic chemical mixtures. 
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There is a recent renewal of interest in stochastic modeling of chemical systems which came with 
the recent realization of importance of noise in cellular information processing. At the level of a 
single cell, number of molecules involved in some processes can be very small and concentrations 
are described as nano-molar J2, 3]. In addition to that, different processes are also characterized by 
significantly different times scales J^]. 

Presence of this time-scale separation and highly different copy numbers of molecular species 
usually complicates the study of biological processes with computer simulations. There is an obvi- 
ous need for computationally tractable stochastic models on a macro-scale that can provide insights 
into joint, qualitative, effects arising from interaction of several sub-networks. In deterministic sys- 
tems of ordinary differential equations, time-scale separation is usually related to the concept of 
stiffness. It is obviously hard to define the same concept in case of the stochastic systems [5]. 

In spite of these obvious complications some progress has been made in modeling of biochemical 
networks which express the separation of time-scales. One difficulty is heterogeneity of simulation 
techniques used for simulation of ODEs/SDEs and stochastic simulation algorithm. One strategy 
exploited in the literature [5, 6] is based on grouping together of reaction events taking place in a 
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single reaction channel in a fast succession and applying diffusion approximation 7]. In [8] Rao 
et.al. discuss a computational approach for performing elimination of the fast species based on rapid 
equilibrium in the limit of the infinite time-scale separation. This method was termed quasi-steady 
state approximation (QSSA). A somewhat similar approach is taken in |9]. Formally, this method 
stems from the classical deterministic QSSA applied to the chemical master equation itself rather 
then to the (stochastic) differential equation underlying the dynamics of the state vector (numbers 
of molecular species). The method developed by Cao et al. in 1 10] can be viewed as generalization 



of approach of Rao et.al. 1 8] but still have the limitations of being derived through the application 
of deterministic techniques and assumptions to the chemical master equation. It also assumes that 
averaging procedure can be done by solving the system of algebraic equations for the expectations 
of the fast variables given slow, termed in 1 1(11 as a virtual fast process. We note here that studies of 
stochastic dynamics of diffusion-type processes evolving on different time scales were pioneered 
by Bogolubov, Khasminski and Freidlin and we refer the reader to monographs 12, 13]. 
This paper has two purposes. First, we present the formulation of stochastic reaction dynamics of 
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reaction network consisting of two subnetworks. Compared to many previous results, where usual 
description of stochastic reaction dynamics follows the approach based on chemical master equation 
(CME), current publication follows the path-sampling approach and represents the dynamics as a 
jump-type stochastic differential equations (SDEs). 

Second purpose is to provide rigorous procedure for the renormalization of the transition rates of 
slow reactions in the presence of fast ones. Following the picture of the stochastic dynamics devel- 
oped in the first part of this paper, we outline the main guidelines for use of stochastic averaging 
principle including error control analysis. Despite of the recent rebirth of interest to the method of 
stochastic averaging in applications to stochastic chemical kinetics, very few examples deal with 
situations when this procedure might break down. We demonstrate here, in a constructive way, 
how to perform the averaging over fast reaction events and how to obtain the effective slow-scale 
transition rate. 

Organization of this paper is as follows. In the next section we discuss the general probabilistic 
framework for stochastic dynamics of reaction networks and introduce a scheme for the partition of 
species and reactions. In Section|ffl]we investigate the consequences of possible time-scales sepa- 
ration and present a procedure based on renormalization of transition rates. We also put emphasis 
on error analysis, outlining main sources of the numerical error on different steps of the procedure. 
Our paper will end with discussion of examples. 



n. NETWORK PARTITIONING 

We begin our discussion with a general set-up, introducing basic concepts and notation. 

Assume that a well mixed, isothermal system has S different molecular species indexed by 
i = 1 ... S and there are R reaction channels, index by r = 1 . . . R, transforming the molecualr 
composition of these species. For the basic notation and examples we direct reader to \\4, 15]. 
State vector of the system can be represented as following: 

(X,Z) (1) 

where fist part of the state vector Xj £ Z+, i = 1 . . . Sx represents main species while the second 
part Z{ S Z+,i = 1. . . S z represents intermidiate species Zi,i = 1 ... Sz- Total number of 
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all types of species: Sx + Sz = S. Vectors , and u^f z are stoichiometric changes of 
components X and Z if reaction event r takes place. We will not make any assumptions about 
actual number of molecular species of each type, i.e. we will not assume low or large copy numbers. 
We assume, however, that there are three subsets of reactions in the system: 

(i) reactions which transform only species X (we denote this subset Tl\), 

(ii) reactions which transform only species Z (subset TZi) 

(iii) "linker" reactions which mix species X and Z (subset 7^3). 

Each reaction channel can be specified by the transition rates a r (a positive function) which 
describes the probability a r dt of reaction event to take place in the interval of time dt. Transition 
rate a r can be further specified as positive functions of X, Z, or,in general, on both components X 
and Z. Based on the definition of subsets 7^1,2,3 we have: 

a r (X), r G TZi (2a) 
a r (Z), r £K 2 (2b) 
Or(X,Z), reTZ 3 (2c) 

We do not assume specific dependence of a r (-) on the state variables X and Z but usually, in 
the framework of mass action kinetics, it is a product of kinetic rate k r and function h r {-) which 
represents the number of reactive configurations available at a given state X, Z I14I1 . 

There exist different methods to characterize the stochastic chemical dynamics. One of the most 
popular approach is to provide an equation for the joint probability density pt(X, Z), which gives 
all information about instantaneous state of the sy_stem at generic moment of time t. Such equation 



is known as chemical master equation (CME) 



er equation (CME) I14L llol and it has be 
QQy. But even if we can obtain 



16] and it has been intensively described and 



utilized in recent literature |8, 91 11011 . But even if we can obtain 1 17] the solution of CME, which 
is usually a very hard problem even for simple chemical networks, this approach still have certain 
limitations, coming from instantaneous description provided by the density pt{-)- 

To describe the stochastic dynamics of the chemical network one can introduce the set of inde- 
pendent point processes N r (t), N r (0) = representing the numbers of reaction events which took 
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place in channels r G 1Z up to time t and use the mass balance relations: 

X t = X(0) + »? N r(t) + V r Z N r (t), (3a) 

Z t = Z(0) + ]T I/fiV r (t) + Y ^r Z Nr(t), (3b) 

where vectors , zv^ and describe the composition change of the system due to the reaction 
event in the channel r. Average number of reaction events in each reaction channel r G ^-1,2,3 
during the small time interval [t, t + 5t) are proportional to the transition rates (|2ji: 

E(N r (t + St) - N r (t)\X t , Z t ) = a r (X t , Z t )5t + 0(5t 2 ) (4) 

Processes N r (t) can be considered as time-changed, unit-rate independent Poisson processes 

n r (i) B: 

N r (t) = U r (f a r {X t ,,Z t ,)dt') (5a) 
Jo 

Thus, the large class of discrete event systems with totally inaccessible event times can be viewed 
as a standard Poisson process with appropriate change of the time scale: 

-t 



[ a r (X t >,Z t >)dt' 
Jo 



(5b) 



The time change generates path-dependent or self-affecting point processes whose dynamics de- 
pend on the information generated by the arrivals of the process (X t , Z t ) . It is important to take 
into account that the stochastic differential equation does not only introduce the probability distri- 
bution for the pair (X, Z) but also generates a measure on the paths, which contains much more 
information. For almost any realization of the set of 1 ... R standard Poisson processes, U r (t, u>), 



parametrized by the element u> of event space 



18] and any deterministic initial condition the 



solution (X(i, co), Z(t, lo)) is a step-wise stochastic process. 

Note also, that dynamics of each component X or Z is non-Markovian if considered separately 
but the dynamics of the pair (X, Z) is Markovian. 

So far we have introduced only the basic notation: quite generic system of SDEs given by © 
outlined in this section have not invoked any assumptions on particular relations between different 
transition rates a r and was totally based on prior information about existence of two groups of 
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species, i.e. X,; and Zj which uniquely identified the partition of the reactions into the subsets 
Hi, IZ2 and 7^.3 . 

In the next section we consider the particular implication of time-scale separation including the 
extensions of the stochastic averaging principle and diffusion approximation. 

m. SEPARATION OF TIME-SCALES AND ELIMINATION OF FAST 
STOCHASTIC VARIABLES. 

In many situations, dynamics of main species X is propagated via large number of fast transitions 
which transform mainly intermediate species Z. One usually desires to construct an approximate, 
time coarse-grained model, which involve only main species. It is important that approximate prob- 
lem describes the dynamics of the system on a large time scale and thus is more advantageous for 
performing simulations without significant sacrifice in accuracy. This section deals with substitu- 
tion of the original problem with approximate one and demonstrates the form convergence of the 
approximation under certain assumptions. 

We assume that at certain region of state space the following assumption can be made about 
transition rates a r (-): 

a r oc O(l) while a r oc 0(e _1 ) (6) 
where separation of the time-scales is introduced via the small parameter e« 1. Problems of this 
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type are challenge for direct application of Stochastic Simulation Algorithm (SSA) 
they will require the time steps of the order 0(e) with a total computational cost of order e _1 . If 
we want to advance through the time interval [0, t], t ~ 0(1) most of the simulation time will be 
spent on simulation of reaction events with the high intensity (X^ r e7^ 2 ar ^ 

O^r 1 )). We would 

like to find an effective transition rates a r (-) for the "linker" reactions (subset 7^3), which describe 
the transition events of the slow reactions "coarse-grained" over the possible events corresponding 
to the reaction events in subset 7^2 • 
It is instructive to consider a simple reaction scheme involving three species Xi, Zi^ similar to 



one considered in 



Zi ^ Z 2 hx 1 (7) 



where rates ki >2 oc e~ l are parametrized by small e and £3 oc O(l). In this case reactions Zi^Z2 
forms the subset 1Z 2 while reaction Z2^Xi corresponds to the subset 1Z 3 and subset 1Z\ is empty,i.e 
TZi = {0}. Then systems of equations for components (X\, Z\, Z 2 ) is the following one: 

Z lt = Z 10 -N 2 (t)+N 3 (t), (8a) 
Z 2t = Z 20 + N 2 (t) - N 3 (t) - N^t), (8b) 
X u = X 1 (0) + N 1 (t) (8c) 

Presence of the scaling factor e _1 in reaction constants fci 2c" 1 allows us to consider family 
of solutions parameterized by e. We expect Z± j2 to follow adiabatically the X\ t . To make that 
apparent, one can apply the functional law of large numbers to the processes N 2j3 (t) in time interval 
[0, t] (see Eqn. d5ali): 

N 2 (t) - N 3 (t)^ (J k 2 Z lv dt' - j k 3 Z 2t 'dA + (9) 

+4= (w 2 ( f k 2 Z w dt') - W 3 ( f k 3 Z 2r dt')\ , e - (10) 
v e V Jo Jo-i / 

where W 2j3 (-) are two independent Wiener processes |7j]. Since parameter e _1 is large, we can 
conclude that difference 

/ k 2 Z\ s ds — I k 3 Z 2s ds 
Jo Jo 

also converges to zero for times t < e/(k 2 + ^3) in the limit of small e, and we can conclude that: 

sup \k 2 Z w - k 3 Z 2v \ (11) 

0<f<t 

This means that variables Z\ t and Z 2t reach a stationary binomial distribution: 

^\Z U Z 2 \X X ) ocQ Zl (l-a) Z2 , (12) 

Z = Z^O) + Z 2 (0) = Z X + Z 2 , a = r 44" ( 13 ) 

k 2 + k 3 

on the time scale t oc 0(e) while sum Z\ t + ^ changes on the much larger time-scale t > 0(1): 

Z lt + Z 2t ^ Z 1 (0) + Z 2 (0)-N 1 (t), (14a) 
Xu^N^t) (14b) 
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By exploiting the separation of time-scales using the stationary distribution 7r e (Zi, Z 2 \Xi) one 
can replace dynamical quantities f{Z\ t , Z2t, X\t) averaged on the time interval [0, t], k 2 +k 3 < ^ 
t < ^- with their conditional averages: 

f(Z lu Z 2u X lt )^- I f(Z w ,Z 2t ,,X w )dt' « (15) 
1 Jo 

« f(X lt ) = f(Zi,Z 2 ,X lt )n e (Z 1 ,Z 2 \X lt ) (16) 

and eliminate fast variables Z12 from the description even though the total number of molecules 
Z\ + Z 2 may be not a large quantity. Thus, taking /(•) to be the "linker" transition rates 
a r (X, Z), r G 7^-3 one obtains averaged transition rates a r (X) which now depend only on the 
slow variable X. Results of the large deviation theory \vi ] demonstrate weak convergence bounds 
of the original problem with small but non-zero e to the solution of the averaged system. But as we 
mentioned it before, one of the goals of this publication is to analyze and extend averaging process 
to the situation when e may be small, but not 'infinitesimally' small. In the next section HTl Al we 
will try to answer this question. 



A. Renormalization of fast fluctuating reaction rates and reduced evolution equations 

Recall that transition rates a r {-) of a jump Markovian process can be used to describe distribu- 
tions of the waiting times of the reaction events via the survival probability of a given state (X, Z) 
has an exponential form S(t) = e~^r=i a <-(X-,z)t an( j describes probability that no reaction event 
take place in any of 1 ... R reaction channels in time interval [0,t] [20]. 

Consider the first jump time of a particular reaction r in the subset of the "linker" reactions, r r> s 
and first jump times of any reaction in the subset of the fast reactions which we will denote t t - :2 . 
Reaction in the group IZ3 have both types of chemical species (X and Z) as their substrates, that 
means that reaction rates in this subset are fluctuate with fast variables Z. If system is originally 
prepared at the state (Xo, Zq) at t = then at any moment of time t > one is interested in finding 
the probabilities of events {r r ^ > t} and {r r> 2 < t}. In other words one has to find an averaged 
survival probabilities: 

S r (t\X) = P({r r , 3 > t}) = ^exp(- a r (X , Zf,)dt')^ , re7J 3 (IV) 
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Average (■ ■ -)z stands for the average over the possible trajectories of the stochastic process 



11]. 



Z x ([0, t]), Z x = Zq at fixed X which depends on X as on parameter 

Probabilities (fTTb can be used to introduce time-dependent transition rates a r (X, t) which effec- 
tively describe the dynamics for reactions in the groups IZ3. Taking the logarithm of the averaged 
survival probabilities ( fTTt we obtain: 

5 r (i|X) = exp(— / dt' d r (X,t')), 
Jo 

a r (i,X) = -J^ln/exp(- J dt' a r (X ,Z t < 

Equations 1 1 8al constitute one of the main results of the paper. In the field of chemical kinetics a 
similar methodolo gy i s known under the label of the "rate dependent processes with dynamical dis- 



(18a) 
(18b) 



order" UA\22, 



23 



24 



25 



2611 where it describes the influence of the non-equilibrium environmen- 



tal degrees of freedom on transport and kinetic properties. Similar approach was used to describe 
juantum dynamics in fluctuating environment [27]. Using the procedure of the cumulant expansion 



quai 
Il6, 



2811 we can obtain the following interrelationship between a r and the multi-point cumulants 



Cr (ti, . . . , t m |X) of the functions a r (X, Z x ), taken at different temporal points t±, . . . , t m : 



5 r (t|X) = exp 



/ dtt... dt m C^ m \t!,...,t m \X) 
^0 ml J * J » 



(19) 



I i\m— 1 rt rt 

.{t,X) = (a r (X,Z t )) z + — ZTi / dt x ... j dt m C { T\t u ...,t m \X) (20) 



m> 2 



m 















20]. 



Renormalized transition rates a r (t, X) provide so-called semi-Markov approximation 
Term "semi-Markov" generally describes non-Markov processes since the statistical properties of 
the waiting times can not be provided only by average rate of the process but all the multi-time 
joint probability distributions for the considered process must be considered. Note that in our case 
effective rate a r depend on the statistics of fluctuations of fast variables Z through the cumulants 

Cr \t\, ■ ■ ■ , t m |X). 

Taking a leading term at e — > 0, which sometimes called Markovian limit, we formally arrive to 
the results of the QSS Approximation 



Or(X,t) = C^ 1) (t|X) = limV a r (X,Z)TT e x (Z) 

e— >0 ^ — J 



(21) 



10 



where average is taken over the invariant measure 7r e (Z|X) of the fast process Z% at fixed X. Note 
that at this level a r does not depend on time and correspond to the single exponential form of the 
survival probability. This level of approximation corresponds to the assumption that at fixed X all 
state space of Z is totally accessible, i.e. ergodic 1 11] and for any function /(•) : U 12 — ► R: 

/"(X) = lim t- 1 / /(X, Z x s )ds = lim V /(X, Z> e (Z*|X) (22) 

There is a general Jensen inequality , which gives the relationship between the mean value of a 
convex function of a random variable an the value of this function when its argument equals the 
mean value of the random variable. According to this inequality: 

S r (t\X) > exp (- j dt'c!p{t'\X)\ (23) 

Application of this inequality leads to the important conclusion that mean field rate d2TT) is larger 
then the rate given by ST% . The exponential and non-exponential structure of the averaged survival 
probability is governed by the hierarchy of the time scales of the dynamics of Z t at different values 
of X. If dynamics of Z is complicated and exhibit metastability at some values of X then Marko- 
vian approximation|^is no longer holds and additional corrections corresponding to the high order 
cumulants must be taken into consideration. Correction to the Markovian approximation based on 
the second order cumulants is: 

Aa r (i,X) = - / dt' <^ 2) (M'|X), (24a) 
Jo 

C®(t, t'\X) = <a r (X, Z?K(X, Z?,)>z " MX, Zf)) z (a r (X, Z%)) z = ((a r (X, Zf)a r (X, Z* v ))) z 

(24b) 

(2) 

The simples assumption for the time dependence of the cumulant Cr is exponential decay: 

Cf)(M'|X) = Kexp{-K(X)\t-t'\) (25) 

where /t(X) _1 is a characteristic relaxation time of the regression of fluctuation of species Z and 
K = ((Aa 2 ,(X, Z)))z- In this case correction to the Markovian term is given by: 

d 

Aa r {t,X) -Kk- 1 ^)— (i-«;- 1 (X)[l-exp(-K(X)i)]) 

Correction to the Markovian approximation given by d24al) is exact for the Gaussian and Markov 
process since the only possible expression for the correlation function of a stationary Markov and 
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Gaussian process is the exponential of a form d25t . It is also interesting to note that correlation 
correction (I24al) generally decreases the transition rate. This is a result which can not be obtained 
using only straightforward averaging method presented in publications Jgl l(3l . 

Note that in general relations (II Sal) can be viewed as a type of fluctuation-dissipation relations; 
they connect the effective dissipation rate in the slow coarse-grained dynamics and statistics of 
fluctuations of the fast reaction events given by the cumulants Cr {t\ , . . . , t m |X) . 



IV. COARSE-GRAINED DYNAMICS AND ERROR CONTROL 

Given the renormalized survival probabilities and transition rates at different points of state space 
of main species X: 

a r (t, X) = a r (X), r 6 IZi 
stochastic dynamics of the main species X can be formulated in the straightforward way, similar to 



the stochastic simulation algorithm (SSA) 
an overall survival probability: 



151 . At the time point t = state Xq we consider 



S(t\X ) = J] ^(*|X Q ) ( 26 ) 

and define a jump moment of the slow process as a first time t\ when 5(t|Xo) crosses the value u, 
where the last one is a random number uniformly distributed on the interval (0, 1) 12911 : 

n = inf{* > 0|S(t|Xo) < u}, «£W(0,1) (27) 

Post-jump transition kernel is defined by the vector of transition probabilities 

o r (ri,Xo) 

q r = ^ —, — T) r e Ki >3 (28) 

Lr'e^, 3 Or'(ri, A ) 

i.e. reaction event r* G 72-1,3 is selected based on the vector q r and current state is updated: 

X Tl = Xo + V r - , t\ = T\ , 

Then the same procedure is performed starting at the state X Tl with generation of the interval T2 
from the survival probability 5(t|X Tl ) and new state X Tl+T2 and so on. As a result one obtains a 
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coarse-grained trajectory: 

n 

(t n ,x tn ), t n = Y^n (29) 

1=1 

Question about the overall accuracy and the error control is a delicate question. Below we de- 
compose the overall error of the method it into the following main factors: 

1. Error in approximating by coarse grained dynamics: 

ei = sup EflXt - X t | 2 ) 

0<t<T 

assuming that transition rates a r {-) can be obtained without error. 

2. Approximation and Monte Carlo error e2 of a T {-) via the finite number of samples represent- 
ing the dynamics of Z t at fixed X. 

Below we discuss step by step leading terms in e\, e-i- 

Estimation of the error e\ is related to the answer on the following question: what possible error 
is introduced while performing averaging of rates of reactions in the subsets 72-1,3 at fixed X? 

It is not hard to see that this error is proportional to the probability of the event that minimal jump 
time over the reactions in group 72 1 U 72-2 is smaller then t while the minimal jump-time of reaction 
in the group 72-3 is larger then t: 

S r (t\Z) = P ({mm r r 9 < t] U { min r r > t}] = /exp(- / dt' a r (X t , , Z ))\ , r G 72 3 
V ren 3 ' re%. 3 J \ Jo I x 

(30) 

where average (. . .)x is taken over trajectories Xf at fixed Z It is not hard to see that this probability 
is exponentially small, i.e. oc exp(— e" 1 ^^) in the limit e — > 0. 

Error e2 depends on the number of cumulants we have included in Eqn. (fT9b and cumulant of 
order m usually gives contribution proportional to e m . In Appendix we outline the exact method for 
calculation of the renormalized survival probability based on eigenvalue decomposition of certain 
linear operator which is a practical approach in situations when state space of the variable Z is not 
very large. 
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v. EXAMPLES 



We now present a simple intuitive example to show that exponential or non-exponential structure 
of the averaged survival probability is governed by the relationship between time-scales of "fast" 
and "slow" species. Assume that for some reaction channel 



X + Z + ... 



(31) 



rate a r (X, Z) = k r h r (X)h' r (Z) jumps reversibly between two values a r (X, 0) and a r (X, 1) with 
the stochastic dynamics of Zt governed by simple master equation: 

^ >l\ V) I 

(32) 



^(0) 



V 



t(l) 



V 




-k 01 kio 
fe i -kio 

Equation (1321 describes the switching transitions between the two states and 1. Assum- 
ing that state of variable Z is prepared according to the equilibrium density it = (7ro,7Ti) = 
( k i+k 10 ' fcof+fcio )' ^ e avera £ e survival probability (e~^o a T {x,z* )dt' ^ can De obtained as follows 
(see also Appendix section for the general computational framework): 





(l\ 


T 


f ( 


S,.(i|X) = 


exp 


t 








\ V 




(33) 



-a r (X,0) - k i ho 

k i -a r (X, 1) - kio. 

This result is very similar in nature to the result obtained in I23I for the case of identical transition 
rates. Remarkable and simple result outlined by Eqn. d33l > allows us to capture in essence regimes 
corresponding to the different ratios of the time-scales: a r <C {kio + ^01) an d a r > (kio + ^oi)- 
First regime (a r -C (kio + ^01)) corresponds to the situation when transitions between different 
states of Z happens much faster then the average rate a r (X, 0), a r (X, 1) of the "linker" process 
and represents the mean-field (MF) regime. In this case dependence of ln(S r (t)) on time t can 
be very well characterized as linear Fig.|2l Not surprisingly, other regime, i.e. a r S> (Alio + koi) 
can be characterized as gated: in this case effective transition rate a r is characterized by the rate of 
switching of Z: koi + kio- 

Figure [5] demonstrates influence of the second order correlation correction Eqn. (l2"5t : 
Aa r (i,X) = 7Ti7To^(l — |(1 — e _K *)), k = koi + kio which fluctuation correction to the ef- 
fective rate a r (-) 



14 



Dependence of survival probability 5 r (i|X) in the example of a two-state system can be shown 
to be non-exponential on the longer time scale but ]n(S r (t)) behaves linearly with time at small 
times t < l/a r (X, •). 



Interesting case of non-exponential relaxation kinetics, and specifically non-exponential kinetics 
at small times can be presented by the following example. Consider a fast reaction given by the 
dimerization reaction: 

S + S S S 2 (34) 

k 2 K eq 

where the fast variable Z t is the number of reaction event which took place up to time t which 
relates the numbers of monomers and dimers with the total number of molecules N m = 2S + S2 
in the following way: 

S = N rn - 2Z, S 2 = Z (35a) 

and a "linker" process is described by the relaxation rate depending on the number of dimers X in 
the following way: 

ar{x > z) = YTx m 

Current value X serves as an activation threshold: at small values of X (X oc 1) only small values 
of Z contribute to the effective rate but probability that Z takes values away from its average are 
exponentially suppressed (Fig.0Jl. On the contrary, if X is large i.e. X « Y2z ~k{Z\X)Z then 
rate given by Eqn. d3oT i depends on the typical value of Z and S r (t\X) manifests time dependence 
similar to the previous example. One can see that this relaxation process shows non-exponential 
time dependence at small times due to the fact that process Z t rarely visits the states contributing 
to the maximum of the relaxation rate given by Eqn. (I36t . We investigate the dependence of the 
averages survival probability on the level of activation threshold X and value of the equilibrium 
constant K eq . Results presented on the Fig.[5]show non-exponential behavior of averaged survival 
probability for the system at small times t. It is evident that non-exponential behavior of S r (t\X) 
is less pronounced for large values of X( X (Z)). 

Eigenvalue-eigenvector decomposition and calculation of expansion coefficients was performed 
via standard routines of LAPACK library available at |http: //www.netiib.orgl . 
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VI. DISCUSSION AND CONCLUSIONS 

Let us summarize the main aspects of this paper. We have studied reduction approach to elimi- 
nate a fast intermediate in the chemical reaction network. To develop this method it is important to 
consider the time coarse-grained transition rates. We have discussed the limitations of the principle 
of stochastic averaging and its possible extensions through the rigorous technique for construction 
of the effective transition rates. We outline the procedure for re-normalization of the transition rates 
and construction of the effective Markov chain for the slow reactions. The merit of the present ap- 
proach is that it is based on a conceptually transparent probabilistic approach involving the waiting- 
time distribution.Technique itself resembles a non-Markovian generalization of the Kubo-Anderson 
theory of stochastic modulation. Our study clearly indicates importance of details of the statistical 
structure of averaging process. 
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vni. FIGURES 



a r (X,0) 



^0 

a r (X,1) 



FIG. 1: Schematic representation of the two-state model. Relaxation rates a r (-) depend on both state Z and X and can 
be quite general. 



19 




FIG. 2: Time dependence of survival probability S r (i) for different ratios of transition rates e — a r (X, l)/(fcoi + fcio) 
for the system with a r (X, 1)^0 and a r (X, 0) = 0. 



APPENDIX A: CALCULATION OF AVERAGED SURVIVAL PROBABILITY 

Calculations of averaged survival probabilities S r (t\X.) requires, in general, the calculation of 
the cumulants of different order m but for some simple cases it can be obtained exactly. This 
is possible for the class of systems which have only finite number of accessible states of the fast 
variables. 
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a r t 



FIG. 3: Time dependence of the survival probability S r (t) calculated with mean-filed (dotted line) approximation and 
second cumulant correction (dashed line) compared to exact dependence (solid line). 

One can study the distribution of values S of the functional 

exp(- [ a r {Z t ,)dt'), (Al) 
J o 

where we have omitted the current state X to simplify the notation. We introducing the joint 



21 



trejactory 



N 




N 



t,time 



probability 



FIG. 4: Trajectory and probability density of the process Z(t). Dotted and dash lines on the probability plot correspond 
to the profile of the relaxation rate a r (X, Z) for different X. 



probability density q(S, Z, t) of the random variables S and Z 

dq(S,Z,t) d 



dt 



a r (Z) — (Sq(S,Z,t))+ 



+ ( a AZ - u r ,))q(S,Z - iv, i) - a r ,(Z))q(S,Z,t)) 







(Z)-(Sq(S, Z,t))+J2 W ZZ 'q(S, Z, t) 



Average survival probability can be expressed following: 



S r (t) = J2 [ Sq(S, Z, t)dS = ^ q r {Z, t) 



z " u z 
and q r (Z, t) is governed by the following master equation: 

dq r (Z,t) 

= —u r \£j )q T \zj : i) t 

Z' 



dt 



T (Z)q r (Z,t) +Y, W zz>q r (Z',t) 



(A2a) 
(A2b) 

(A2c) 



(A2d) 



(A2e) 



One can find an averaged survival probability via eigenvalue-eigenvector {A, V\(Z)} decomposi- 
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FIG. 5: Time dependence of the survival probability S r (t\X) for the system where dimerization dynamics of the fast 
variable Z is described by parameters N m — 200, fci = 1.0, ki = 10.0, K eg = 10 2 . Plots are shown for values 
of X — 1 and 50 clearly manifest non-exponential character of the relaxation process at small time for low values of 
X. Note that kinetics is non-exponential on time larger then characteristic scale t non - exp « 0.02 of fluctuation of Z ( 
k± 1 (N/2) 2 w 10~ 3 ) i.e. on the relevant for time-coarsening interval. 
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tion of the linear operator W zz' — a r {Z)5zz''- 

Z A 

where coefficients c\ correspond to the decomposition of the invariant probability n(Z\ •): 

vr(Z) = ^c a Va(Z) (A4) 

A 



